Efficient tunable generic model for fluid bilayer membranes 
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We present a model for the efficient simulation of generic bilayer membranes. Individual lipids are repre- 
sented by one head- and two tail-beads. By means of simple pair potentials these robustly self-assemble to a 
fluid bilayer state over a wide range of parameters, without the need for an explicit solvent. The model shows the 
expected elastic behavior on large length scales, and its physical properties (e. g. fluidity or bending stiffness) 
can be widely tuned via a single parameter. In particular, bending rigidities in the experimentally relevant range 
are obtained, at least within 3 — 30fcBT. The model is naturally suited to study many physical topics, including 
self-assembly, fusion, bilayer melting, lipid mixtures, rafts, and protein-bilayer interactions. 

PACS numbers: 61.20.Ja, 81.16.Dn, 82.70.Uv 



Lipid molecules in aqueous solution spontaneously assem- 
ble into bilayer membranes. In biological systems, such mem- 
branes are involved in tasks over an extraordinary range of 
length scales, from transport of water and ions at the scale 
of nm, up to phagocytosis, amoebal motion and cell bud- 
ding at the scale of jum |]]|. Computer simulations designed 
to understand some aspects of this structural and functional 
range must accordingly be tailored to the specific length and 
timescales involved. Techniques that probe both the small- 
est J3l and largest [33 of these length and time scales are now 
comparatively well established; however, accessing interme- 
diate regimes has proven far more difficult, and it is only re- 
cently that significant progress has been made in this regard. 
The need for a comprehensive suite of techniques to study 
lipid bilayers at the mesoscale is highlighted by the sheer 
number of relevant problems in this regime, which include 
viral budding, raft formation, fusion, phase separation of mul- 
ticomponent systems and protein sorting during vesiculation. 

Most existing approaches to mesoscale simulation employ 
coarse grained lipids and require explicit solvent particles to 
stabilize the bilayer. This strategy is convenient and natural, 
yet it comes at a high price: Already for small flat systems 
the solvent accounts for most of the computational effort, but 
the problem gets significantly worse when three dimensional 
objects such as vesicles are to be simulated. Membrane and 
solvent play the role of surface and bulk, respectively, hence 
the solvent degrees of freedom vastly outnumber the lipids 
even for rather modest sized vesicles. An obvious solution 
to this problem is to replace the solvent by effective lipid in- 
teractions. Given the great success of this approach in poly- 
mer physics it is perhaps surprising that solvent-free bilayer 
simulations have so far failed to find widespread acceptance. 
The problem appears to be that naive choices for interparticle 
potentials (e. g. Lennard- Jones) do not lead to a fluid bilayer 
phase but only to ordered "solid" bilayers at low temperature 
and low density phases at high temperature. Many attempts 
to obtain a broadly stable fluid phase have been made with 
varying degrees of success. So far, however, none of these has 
resulted in a model that is sufficiently robust, simple, or ver- 
satile for general use. For example, the original solvent-free 
model of Drouffe et al. 01 and later modifications by Noguchi 



0] rely on density dependent interaction potentials to stabi- 
lize the fluid phase. But the multibody nature of interactions 
poses serious problems for interpretation and measurement of 
thermodynamic quantities. Other models do not exhibit the 
crucial property of unassisted self-assembly |6, 7] and require 
the use of angular dependent potentials |7] or a large set of 
highly tuned interaction parameters 0], Thus, there remains 
a clear need for an efficient, robust and tuneable solvent-free 
bilayer model. 

In this letter we present a very simple method for simu- 
lating solvent-free fluid bilayer membranes that is based on 
pair potentials, is highly robust and tuneable, and reproduces 
macroscopic properties of real bilayers. The key ingredient is 
an attractive potential between lipid tails for which the range 
of attraction, w c , can be varied. First we map out the prop- 
erties of a lipid system as a function of w c and temperature 
and find that a sufficiently large w c is the key to obtaining a 
fluid phase. Moreover, w c can be used to tune bilayer prop- 
erties and to construct multi-component systems with inter- 
facial tension. This is demonstrated in the final part of our 
Letter where we examine the kinetics of domain formation in 
a two-component vesicle. 

Let us now describe a model that demonstrates the principle 
mentioned above. Each lipid molecule is represented by one 
"head" bead followed by two "tail" beads. Their size is fixed 
via a Weeks-Chandler- Andersen potential 



KepM) = 



, r > r c 



(1) 



with r c = 2 1 / 6 fr. We use e as our unit of energy. In or- 
der to ensure an effective cylindrical lipid shape we choose 
&head,head = fyiead.taii = 0.95 a and fc» t aii,taii = cr, where o~ is the 
unit of length. The three beads are linked by two FENE bonds 



VbondM = -^bond^L l°g [l ~ (r/roo) 2 ] 



(2) 



with stiffness fcbond = 30 e/ a 2 and divergence length r M = 
1.5 a. Lipids are straightened by a harmonic spring with rest 
length 4cr between head-bead and second tail-bead 



HendM = ^k bead (r - 4cr) 5 



(3) 
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FIG. 1 : Phase diagram in the plane of potential width w c and temper- 
ature at zero lateral tension. Each symbol corresponds to one simu- 
lation and identifies different bilayer phases: x: gel; •: fluid, +: 
unstable. The dashed lines are merely guides to the eye. The inset 
shows the pair-potential between tail lipids (solid line) and the purely 
repulsive head-head and head-tail interaction (dashed line). 

which corresponds in lowest order to a harmonic bending po- 
tential ^fcbendC 2 $ 2 for the angle n — i? between the three 
beads. We fixed the bending stiffness at Amende 2 = 10 e. Fi- 
nally, all tail beads attract each other according to 

{-e , r < r c 

-e cos 2 , r c <r<r c + w c , (4) 

, r > r c + w c 

This describes an attractive potential with a depth of e which 
for r > r c smoothly tapers to zero. Its decay range w c is the 
key tuning parameter in our model. 

We performed Molecular Dynamics (MD) simulations us- 
ing the ESPResSo package |8]. A fixed number of lipids was 
simulated in a cubic box of side length L subject to peri- 
odic boundary conditions. The canonical state was reached 
by means of a Langevin thermostat |9] (with a time step 
St = 0.01 t and a friction constant T = t _1 in Lennard- 
Jones units). If needed, constant tension conditions were also 
implemented via a modified Andersen barostat 1 10] (with a 
box friction Tbox = 2 x 10 _4 t _1 and box mass within the 
range Q = 10~ 5 . . . x 10~ 4 ). 

One of our principal objectives is to map out the conditions 
under which the fluid bilayer is stable. We identified the fluid 
phase in two different ways. First, a box-spanning bilayer 
was pre-assembled from 1000 lipids, and its equilibration un- 
der zero lateral tension was attempted (requiring - if success- 
ful - box lengths of L « 25 a). Three qualitatively different 
outcomes were observed: (i) At sufficiently low temperature 
the bilayer adopted a "gel" phase; (ii) within a more elevated 
temperature range a fluid phase can be reached, provided 
w c > 0.8 a; (Hi) at sufficiently high temperature a bilayer 
under zero tension always fell apart. Fluid and gel phases 
are already clearly distinct under visual inspection (long range 
order in the gel phase and none in the fluid). Across the gel- 



fluid boundary flUl we observed a sharp increase in in-plane 
diffusion constant D, an abrupt decrease in orientational order 
5 = i(3(a i -n) 2 — 1}^ (where is the vector along the axis of 
the I th lipid and n is the average bilayer normal) and the emer- 
gence of a nonzero flip-flop rate rf (the probability per unit 
time that a single lipid changes its monolayer). Typical values 
for these parameters along the k^T = Lie isotherm are D = 
0.06 - 0.03cr 2 /r, S = 0.5- 0.8 and r f = 2 - 90 x IO^t' 1 . 
For uu c = 1.6cr and k^T = Lie we also measured the bilayer 
rupture tension £ r 4mN/m and the compressibility mod- 
ulus at zero tension, JC w 50mN/m (mapping to real length 
scales by assuming a bilayer thickness of 5nm). These quan- 
tities are within their respective measured ranges for synthetic 
and biological membranes, £ = 0.1 — 12mN/m II 211 and 
K, = 50 — 1700mN/m (depending on cholesterol content) 
1131 . Second, lipids were simulated at constant volume, but 
starting from a random "gas" configuration. Under all con- 
ditions which previously gave stable tensionless membranes, 
a bilayer patch quickly self-assembled which at the right L 
could zip up with its open ends to span the box. If the box 
was too big, the patch either remained free, or (sometimes) 
closed upon itself to form a vesicle. Box spanning bilayers 
could also occur somewhat above the evaporation boundary 
of Fig. Hi 14], indicating that this line should be viewed as the 
location where the rupture tension for a bilayer approaches 
zero. We finally remark that upon approaching this line from 
below, the bilayer becomes increasingly disordered, and flip- 
flop rates and diffusion constants increase strongly. 

Looking at figure ^ we see that the temperature range over 
which a fluid bilayer is stable increases as w c is increased and 
that for relatively narrow potentials w c < 0.7a the fluid region 
disappears completely 1 15]. It is noteworthy that these general 
features are not restricted to the present functional form of 
Eqn. @. Indeed, we have also obtained a qualitatively similar 
phase diagram for a Lennard-Jones like potential with variable 
width. It should therefore be emphasized that it is not the 
precise functional form of our tail attractions that is important 
for the stabilization of fluid bilayers but rather the length scale 
over which these attractions are effective. 

Next we illustrate that the fluid bilayers approach the cor- 
rect elastic continuum limit. If one expands the bilayer in 
modes h(r) = ^ q h q e lq r with q = ^(n x ,nj,), (linearized) 
Helfrich theory predicts the mode spectrum [16] 

where k is the bending modulus and a the lateral tension. Be- 
low the crossover wave vector q c = \J a j k one has {\h^\) ~ 
q~ 2 (tension regime), while (|/i 2 1) ~ g~ 4 holds above (bend- 
ing regime). Once 1/q approaches length scales comparable 
to the bilayer thickness, continuum theory breaks down and 
further effects (e. g. protrusion modes 11711 ) set in. In order to 
see the characteristic q~ A scaling of the bending regime over a 
sufficiently wide range requires q c to be as small as possible, 
hence we simulated at zero tension. Furthermore, to get away 
from microscopic lengths (such as the bilayer thickness) we 



FIG. 2: Asymptotic q~ 4 scaling of the power spectrum (\hq\) for 
the bilayer system with w c = 1.6a and feT = Lie. The inset 
shows the thus measured bending stiffness values at fixed tempera- 
ture (feT = Lie) as a function of potential range ui c - 



took systems four times as big as the ones we used for map- 
ping the phase diagram (4000 lipids, L ~ 50 a). Note that 
reaching the continuum limit in MD simulations is not trivial, 
since the relaxation time of bending modes scales as q~ 4 . 

After setting up the bilayer, we first waited until tension, 
box length, and energy had equilibrated (which took typically 
10 5 t for fluid systems). Then on the order of 100 configu- 
rations separated by 2000t were used to measure the mode 
spectrum. The bilayer mid-plane was identified by tracking 
the tail-beads and interpolating their vertical position onto 
a 16 x 16 grid. Possible stray lipids had to be excluded 
from this procedure. An FFT then yields the power spectrum 
Fig.|2]illustrates (for the system with w c = 1.6ct and 
k^T = Lie) the q~ 4 scaling. Notice that length scales ex- 
ceeding L w 20 cr (z. <?., about four times the bilayer thickness) 
are required to reach the asymptotic regime, outside of which 
a fit to Eqn. (0 and the subsequent extraction of a bending 
constant is meaningless. The inset shows the corresponding 
values of n for the six simulated stable bilayer systems at this 
particular value of the temperature (see Fig. We empha- 
size that its value is easily tuneable over the experimentally 
interesting range, at least from k ss 3 k-gT up to n w 30 k-sT. 

So far we have used the parameter w c only to tune the 
bilayer stiffness. Let us now illustrate another application, 
namely, the possibility of creating mixed lipid systems. We 
consider a simple example in which two lipid types A and B 
are present and where w AA = . A line tension can now be 
generated by choosing the cross interaction w AB smaller than 
the homogeneous ones. Provided a sufficiently small cross 
term has been chosen, A and B domains will form. Figure 
13 illustrates qualitative features of this process for a critical 
quench (A:B = 1:1). First we note that matching domains al- 
ways appear on inner and outer leaflets. As is typical for a 
critical system, domains are initially elongated in shape and 



FIG. 3: Phase separation and budding sequence for a mixed vesicle 
pre-assembled from 8000 A-lipids and 8000 B-lipids at k B T = Lie. 
After an equilibration time of 2000 r, during which w£ A = w^ B = 
Wc B = 1.5(7, the cross term was reduced to w^ B — 1.3(7. Times 
indicated are measured from this point on. 

then coarsen to form circular patches or stripes. In this case 
the domain size and A-B line tension are sufficient to induce 
budding 1 18]. This process has recently been studied in detail 
via DPD simulations 1 19]. For matching lipid proportions in 
outer and inner leaflets, budding was induced by cleavage of 
the domain boundary. We also observed this "flaking mecha- 
nism" for very high line tensions (w AB <C w AA ' BB ), while for 
less pronounced line tensions a more continuous type budding 
occured (see Fig. O, provided the vesicle was big enough. 

Next we investigate the kinetics of domain formation on 
a larger vesicle in the off critical regime. Although the ki- 
netics of domain formation has been extensively studied for 
supported membranes, only a small number of recent sim- 
ulations fl9l Eol EUl and experiments I22I1 have tackled the 
problem for free vesicles. The importance of these studies is 
underscored by the fact that many features of a vesicular sys- 
tem are absent in a flat supported membrane. These include 
bending fluctuations, volume constraints, and kinetic effects 
of curvature-composition coupling 1 23 ] . 

To study the kinetics of phase separation we performed sim- 
ulations using pre-equilibrated vesicles with an A:B ratio of 
3:7 and measured the number n(t) of clusters of A lipids as 
a function of time. This quantity displays two distinct kinetic 
regimes (see Fig.|^l. In the first regime n(t) decays exponen- 
tially, corresponding to a conversion of an initially exponential 
cluster size distribution to an arrangement where meso-sized 
clusters begin to dominate. At later times such clusters display 
an n(t) ~ t e power law decay with 8 w —0.4, which agrees 
with the value expected for coarsening via patch-coalescence 
in the non-hydrodynamical regime 1 24] . This scenario is con- 
firmed visually, as well as by characteristic jumps in the (size 
weighted) average cluster size. Interestingly, we do not ob- 
serve evidence for a Lifshitz-Slyozov type ripening, for which 
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FIG. 4: Coarsening kinetics for a vesicle composed of 4800 A-lipids 
and 11200 B4ipids at k B T = Lie with self interactions ro c M = 
iuj? B = 1.7a. A cross term of w^ B — 1.5a was imposed at t = Or. 
The number n(t) of clusters of size three or greater was calculated 
as an average over five independent runs. After an initial exponential 
decrease, a power law with exponent 8 ~ —0.4 is observed. 

one would have 9 = —2/3 |25|. 

Having demonstrated the key features of the present model, 
including its tunability and its application to multi-component 
systems, we now turn to performance aspects. Although one 
of the clear advantages of our approach is speed, it is difficult 
to make meaningful comparisons across different computer 
architectures and implementations. One basic quantity that 
provides a rough implementation-independent comparison is 
the particle number. Simulating a vesicle with 16000 coarse 
grained lipids, as in our example, would require roughly 25 
times as many solvent particles [ 19, 20]. (Actually this factor 
scales with the vesicle radius). Obviously, one must also ac- 
count for the time-step which is often chosen somewhat larger 
in DPD simulations. Still, this leads us to an at least 5-fold 
speed-up for the present method over DPD, and a crude com- 
parison based on direct CPU time usage agrees with this esti- 
mate lErill . 

In summary, we have presented a method that is fundamen- 
tally different from existing coarse grained lipid membrane 
simulations that use explicit solvent. Our method should be 
seen as complementary to these techniques since it presents a 
compelling speed advantage, especially for systems in three 
dimensions (vesicles, bicontinuous phases). It does not nat- 
urally include volume constraints or hydrodynamics. Indeed, 
the model is inspired by the great success of simple "bead- 
spring" models used to study polymer systems in the Rouse 
regime. We have presented Langevin Dynamics simulations 
here, but the simplicity of the model and concept easily per- 
mit implementation of other integrators or even Monte Carlo. 
Volume constraints should be relatively simple and efficient 
to include via the concept of a "phantom solvent" fzfa . The 
particle based nature of the model and its tunable interactions 
allow for great flexibility such that components with different 



pre-determined stiffness and lipid shape present no difficulty, 
and topological changes occur naturally. Despite its astound- 
ing simplicity our approach allows one to capture all these es- 
sential aspects of lipid bilayer physics within a single unified 
framework. 

We thank Oded Farago, Friederike Schmid, Hiroshi Nogu- 
chi, and Bernward Mann for valuable discussions. 



[1] H. Lodish, A. Berk, S. L. Zipursky, P. Matsudaira, D. Balti- 
more, J. Darnell, Molecular Cell Biology (Freeman & Com- 
pany, New York, 2000). 

[2] D. P. Tieleman, S. J. Marrink, and H. J. C. Berendsen, Biochim. 
Biophys. Acta 1331, 235 (1997); S. E. Feller, Curr. Opin. Col- 
loid Interface Sci. 5, 217 (2000); L. Saiz, S. Bandyopadhyay, 
and M. L. Klein, Biosci. Rep. 22, 151 (2002). 

[3] G. Gompper and D.M. Kroll, Triangulated-Surface Models of 
Fluctuating Membranes, in: Statistical Mechanics of Mem- 
branes and Surfaces, 2 nd ed., ed. by D.R. Nelson etal., (World 
Scientific, Singapore, 2004). 

[4] J.-M. Drouffe, A. C. Maggs, and S. Leibler, Science 254, 1353 
(1991). 

[5] H. Noguchi and M. Takasu, Phys. Rev. E 64, 041913 (2001); 

J. Chem. Phys. 115, 9547 (2001); Biophys. J. 83, 299 (2002); 

Phys. Rev. E 65, 051907 (2002). 
[6] O. Farago, J. Chem. Phys, 119, 596 (2003). 
[7] G. Brannigan and F. L. H. Brown, J. Chem. Phys. 120, 1059 

(2004); G. Brannigan, A. C. Tamboli and F. L. H. Brown, J. 

Chem. Phys. 121, 3259 (2004). 
[8] http://www.espresso.mpg.de 

[9] G. S. Grest and K. Kremer, Phys. Rev. A 33, 3628 (1986). 
[10] A. Kolb and B. Diinweg, J. Chem. Phys. Ill, 4453 (1999). 
[11] In the interest of brevity we have omitted a detailed analysis of 

the gel-liquid transition and also of possible subdivisions within 

the gel phase itself. These issues will be treated elsewhere. 
[12] C. E. Morris and U. Homann, J. Membrane Biol 179, 79 (2001). 
[13] D. Needham and R. S. Nunn, Biophys. J. 58, 997 (1990). 
[14] This is most likely a finite-size effect (Brannigan, Tamboli, and 

Brown |7] observe a similar ensemble dependence). 
[15] The LI potential corresponds to w c ~ 0.7 a, and has indeed 

historically proven to be unsuccessful in creating fluid bilayers. 
[16] U. Seifert, Adv. Phys. 46, 13 (1997). 

[17] R. Lipowsky and S. Grotehans, Europhys. Lett. 23, 599 (1993). 
[18] R. Lipowsky, Biophys. J. 64, 1133 (1993). 
[19] S. Yamamoto and S. Hyodo, J. Chem. Phys. 118, 7937 (2003). 
[20] M. Laradji and P. B. Sunil Kumar, Phys. Rev. Lett. 93, 198105 
(2004). 

[21] P. B. Sunil Kumar, G. Gompper and R. Lipowsky, Phys. Rev. 

Lett. 86, 3911 (2001). 
[22] S. L. Veatch and S. L. Keller, Phys. Rev. Lett. 89, 268101 

(2002); Biophys. J. 85, 3074 (2003); T. Baumgart, S. T. Hess 

and W. W. Webb, Nature 425, 821 (2003). 
[23] T. Taniguchi, Phys. Rev. Lett. 76, 4444 (1996). 
[24] K. Binder and D. Stauffer, Phys. Rev. Lett. 33, 1006 (1974). 
[25] D. A. Huse, Phys. Rev. B. 34, 7845 (1986). 
[26] Our 16000 lipid vesicle system takes ~ 58 cpu hrs to complete 

2000r on a 2.0GHz AMD Opteron246. This compares with ~ 

300 cpu hrs for a recent DPD simulation |.2y]. 
[27] O. Lenz and F. Schmid, J. Mol. Liquids 117, 147 (2005). 



